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Abstract 

We  study  nontrivial  applications  of  the  reduced  basis  method  (RBM)  for  elec¬ 
tromagnetic  applications  with  an  emphasis  on  scattering  and  the  estimation 
of  radar  cross  section  (RCS).  The  method  and  several  extensions  are  explained 
with  two  examples  with  different  characteristics.  Parameters  that  are  allowed  to 
vary  within  the  model  include  frequency,  incident  angle  and  measurement  angle 
as  well  as  the  geometry  of  the  scatterers.  With  appropriate  applications  of  the 
empirical  interpolation  method  (EIM),  transformation  of  the  domain,  configu¬ 
ration  of  perfectly  matched  layer,  exponential  convergence  of  the  reduced  basis 
solution  over  the  entire  parameter  domain  is  achieved.  Moreover,  we  demon¬ 
strate  that  this  approach  allows  for  the  effective  capture  of  the  critical  behavior, 
in  this  case  through  shapes  that  minimize  scattering.  This  further  highlights 
the  robustness  and  quality  of  the  greedy  approximation  and  the  reduced  basis 
method  approach. 

Keywords:  Reduced  Basis  Method,  Electromagnetic  Scattering,  Radar  Cross 
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1.  Introduction 

In  many  applications  such  as  computational  optimization,  control  and  de¬ 
sign,  one  is  required  to  rapidly  yet  accurately  predict  some  quantities  of  inter¬ 
est  under  the  variation  of  a  set  of  parameters  u  G  V  C  Here,  s^{u)  := 
z/)  G  C  is  an  output  of  interest  with  being  a  functional  and  is  the 
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solution  of  a  parametrized  partial  differential  equation  (PDE).  Let  us  denote 
the  PDE  in  its  weak  form 

a^K(z/),u;z/)  =  r(u;z/),  Vu  G  (1) 

where  a®  and  are  bilinear  and  linear  forms,  respectively,  associated  with  the 
PDE.  We  shall  denote  as  the  space  of  the  exact  solution  u^. 

In  practice,  one  could  use  a  finite  element  (EE)  discretization  to  approximate 
the  solution  {y)  u^{u).  Hence,  given  u  £  V  C  find  {u)  G 

satisfying 

a'^(?x'^(z/),  u;  z/)  = /■^(u;  z/),  Vu  G  (2) 

An  approximate  value  of  the  output  of  interest  can  be  thus  computed  by  {y)  := 
0 {u^ {y)]  u)  G  C.  Here  X-^  is  the  finite  element  space  approximating  X^  with 
dim(X-^)=  A/",  •;  •)  and  m-^ {•]  •),  m  G  {/, /}  are  computable  approxima¬ 

tions  of  a(',  •;  •)  and  •),  m  G  {/,  /},  respectively.  We  assume  {u)  provides 
a  reference  solution,  referred  to  as  the  truth  approximation,  that  is  accurate 
enough  for  all  u  G  V.  To  ensure  this,  one  must  usually  choose  a  very  large  Af. 
As  a  result,  the  cost  to  solve  for  the  truth  approximation  is  likely  very  high, 
especially  when  the  solution  is  needed  for  many  instances  of  z/,  in  which  case 
the  cost  becomes  prohibitive. 

The  RBM,  introduced  in  [15,  9],  offers  a  particularly  well  suited  solution  to 
the  challenges  of  this  “many-query”  scenario,  aiming  to  improve  the  simulation 
efficiency  and  reduce  the  overall  computational  cost.  A  fundamental  observa¬ 
tion  and  assumption  utilized  by  RBM  is  that  the  parameter  dependent  solution 
u^{u)  is  very  often  not  simply  an  arbitrary  member  of  the  infinite-dimensional 
space  X^  associated  with  the  PDE,  but  rather  it  evolves  on  a  lower-dimensional 
manifold  induced  by  the  parametric  dependence.  Therefore,  one  can  expect 
that  as  z/  (g  P  C  varies,  the  set  of  all  solutions  iz®(z/)  can  be  well  ap¬ 
proximated  by  a  finite  and  low  dimensional  vector  space.  In  particular,  the 
RBM  method  proposes  to  approximate  the  solution  of  the  PDE  for  an  arbi¬ 
trary  value  of  the  parameter  z/  G  P  as  a  linear  combination  of  solutions  of 
the  same  PDE  for  other  adequately  chosen  parameters  z/^,  i  G  {1, . . . ,  X},  i.e., 
u^{iy)  =  Ci{v)u^ {vi).  The  truth  approximation  problem  is  then  replaced 

by  the  following  problem:  find  {v)  G  X^  such  that 

a^{u^{v),v;v)  =  y{v;v),  \/v  €  ,  (3) 

where  X^  =  span  {iz-^(z/^),  i  G  {!,...,  X}}.  As  can  be  expected,  the  choice  of 
the  initial  set  of  parameters  used  to  compute  the  basis  functions  appropriately 
is  crucial  for  the  method  to  succeed.  This  is  guided  by  the  combination  of  a 
rigorous  a  posteriori  error  estimators,  also  used  to  certify  the  quality  of  the 
approximations,  and  a  judicious  greedy  algorithm.  We  refer  to  [20,  17,  10,  6] 
and  also  the  review  paper  [19]  and  the  extensive  reference  therein  for  detail. 
Theoretically,  exponential  a  priori  convergence  result  of  the  reduced  basis  ap¬ 
proximation  is  confirmed  for  a  one  dimensional  parametric  problem  [14].  More 
recently,  exponential  convergence  of  the  greedy  algorithm  for  continuous  and 
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coercive  problems  with  parameters  in  any  dimension  has  been  established  in  [5] , 
and  improved  in  [4]. 

The  RBM  becomes  particularly  valuable  when  the  forms  involved  in  the 
PDE  satisfy  the  so-called  affine  assumption,  namely, 

Qa  Qm 

a'^{u,v]v)  =  {v),  (4) 

9=1  9=1 

where  m  G  {/,/}•  In  this  case  the  computational  strategy  centers  around  split¬ 
ting  into  offline  and  online  parts.  In  the  offline  part  (computations  that  are 
performed  once),  the  reduced  basis  space  and  various  elements  needed  for  the 
error  estimation  is  precomputed,  at  a  cost  scaling  with  M.  In  the  online  part, 
where  the  reduced  basis  approximation  is  used  to  approximate  the  solution  for 
any  new  parameter,  the  computational  cost  depends  only  on  N  (the  dimension 
of  the  reduced  basis  space)  and  Qm?  ^  ^  so-called  complexity 

of  the  forms)  but  not  on  M  (the  dimension  of  the  FE  space  where  the  truth 
approximation  is  computed).  Since  N  M  this  introduces  a  potential  for  a 
substantial  computational  saving. 

In  this  paper,  we  consider  problems  of  electromagnetic  scattering  in  which  a 
valuable  output  of  interest  is  the  radar  cross  section  (RCS)  [12],  measuring  the 
reflective  strength  of  a  target  when  illuminated  by  an  electromagnetic  source. 
The  design  and  optimization  of  this  radar  signature  is  naturally  of  great  practical 
and  tactical  importance,  involving  both  material  and  geometric  design  options. 
Other  parameters  to  include  are  incident  and  observation  wave  angles  and  the 
frequency.  In  general  the  linear  and  bilinear  functionals  associated  with  class  of 
problems  do  not  satisfy  the  affine  assumption  (4).  In  consequence,  an  additional 
approximation  based  on  the  EIM  will  be  developed  and  applied  [2,  10].  Reduced 
basis  method  has  been  applied  to  3D  electromagnetic  scattering  problem.  In¬ 
spired  by  the  subdomain  residual  method,  the  authors  of  [16]  proposed  a  faster, 
but  non-rigorous,  error  estimate  to  deal  with  affine  expansions  with  very  large 
number  of  terms. 

What  remains  of  this  paper  is  structured  as  follows.  In  Section  2,  we  state 
the  general  problem  setting  before  introducing  two  specific  scattering  problems 
that  we  shall  use  as  benchmarks  to  illustrate  the  capabilities  of  the  reduced 
basis  method  for  complex  scattering  and  radar  cross  section  prediction  applica¬ 
tions.  These  examples  will  involve  parametric  variations  of  frequency,  incident 
and  observation  angles  and  geometric  variations,  or  combinations  of  these.  In 
Section  3,  we  briefly  explain  the  computational  strategy  and  offer  some  details 
on  the  truth  approximation  solver  and  the  greedy  algorithm  used  to  determine 
the  parameters  enabling  the  construction  of  the  reduced  basis  space  [17,  20,  6]. 
The  two  test  problems  do  not  satisfy  the  affine  assumption  and  we  develop  an 
empirical  interpolation  method  (EIM)  that  facilitates  a  reduced  problem  satisfy¬ 
ing  (4),  ensuring  an  efficient  off-line  on-line  computational  strategy.  Numerical 
results  illustrating  the  superior  performance  of  reduced  basis  method  for  the 
two  non-trivial  examples  are  presented  in  Section  4.  Finally,  we  offer  some 
concluding  remarks  in  Section  5. 
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2.  Presentation  of  the  problem 

In  this  section,  we  describe  the  setup  and  formulation  of  two  non-trivial  bench¬ 
mark  cases  when  the  radar  cross  section  is  sought  as  the  output  of  interest. 
For  the  first  case,  the  angular  frequency  and  angle  of  the  incident  wave  are 
the  parameters,  while  in  the  more  challenging  second  case  the  parameter  is  the 
shape  of  the  scatterer.  As  we  shall  see,  this  necessitates  the  nontrivial  geomet¬ 
ric  transformation  and  the  corresponding  application  of  perfectly  matched  layer 
(PML). 

In  both  cases,  the  electromagnetic  waves  are  TM-polarized,  i.e.,  the  electric 
and  magnetic  fields  satisfy  E  =  (0,  0,  Ez)  and  H  =  0)  in  the  Maxwell’s 

equation.  As  a  source,  we  consider  an  incident  plane  wave 

Eincix)  = 

Here,  uj  is  the  angular  frequency  and  d  =  d{0i)  :=  (— cos(6>i),  —  sin(6>i))  is  the 
direction  of  propagation  that  depends  on  the  angle  of  the  incident  wave. 

The  scattered  field,  generated  by  the  incident  wave,  satisfies  a  set  of  equa¬ 
tions  (we  assume  6  =  /r  =  1  for  the  sake  of  simplicity): 


II 

dHy 

dx 

OH, 

dy 

in  ^4, 

(5a) 

iujHx  = 

_dE, 

dy 

in  f]. 

(5b) 

iUjHy  = 

dE, 

dx 

in  Q, 

(5c) 

Einc 

on  F. 

(5d) 

The  system  is  completed  by  requiring  the  scattered  field  to  satisfy  the  so-called 
Silver-Miiller  radiation  condition  at  the  infinity.  Here,  Q  is  the  domain  outside 
of  the  scatterer  and  T  is  the  boundary  of  the  scatterer.  To  enable  tangible 
computations,  an  exterior  boundary  sufficiently  far  away  from  the  scatterer, 
T far^  is  enforced  to  reduce  the  infinite  domain  to  a  finite  one  and  artificial 
boundary  conditions  will  be  applied  on  the  exterior  boundary  T In  this  work, 
either  the  fields  {Ez^Hx^Hy)  are  assumed  to  satisfy  the  first  order  absorbing 
boundary  condition  (on  the  first  problem)  as 

{Hxn)‘Z  =  ^Ez,  onTfar,  (6) 

or  we  apply  PML  [3,  7,  1]  on  an  annulus  just  inside  of  F far  (on  the  second 
problem). 

Given  the  scattered  field  H  and  E,  we  can  calculate  the  radar  cross  section 
associated  with  this  scattering  process.  We  define  the  surface  currents 

J  =  n  X  H,  and  M  =  —n  x  E, 
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and  let 


i—  n 

F{UJ,  Oi,  Or)  =  i  {-ujfiz  ’J-ujzxM’t)  dC.  (7) 

vSttcj  J 

Here  f ,  the  unit  vector  for  the  observation  direction,  is  equal  to  (cos(6>^),  sin(6>^)) 
with  Or  being  the  observation  angle,  f  is  the  vector  from  scatterer  to  the 
integration  contour  which  can  be  any  closed  contour  surrounding  the  scatterer. 
With  the  linear  output  functional  (7)  computed,  the  RCS  (also  called  bi-static 
RCS)  is  defined  as 

RCS{UJ,  Ou  Or)  =  (8) 

A  case  of  particular  interest  in  practice  is  the  monostatic  RCS  defined  by 
RCS{uj^0i^0r  =  Oi)  .  Often,  these  magnitudes  (quotients  of  powers)  are  mea¬ 
sured  in  decibels  (dB). 

2.1.  Open  Cavity  Benchmark 

In  the  first  problem  to  be  considered,  the  scatterer  is  assumed  to  be  a  perfect 
electric  conductor  (PEC)  open  cavity  given  by  (see  Figure  1  left) 

r  :=  [-1, 1]  X  {-1}  U  [-1, 1]  X  {1}  U  {-1}  X  [-1, 1], 

The  unbounded  nature  of  the  computational  domain  is  enabled  by  a  first  order 
absorbing  boundary  condition  (often  called  the  Silver-Muller  boundary  condi¬ 
tion)  at  the  boundary  T far  =  {x  G  such  that  \x\  =  3}.  The  output  of 
interest  is  the  radar  cross  section  given  by  (8).  For  this  application  we  consider 
three  parameters  {uj^Oi^Or)  £  V  =  [7r/2,57r/2]  x  [0,  27r]  x  [0,27r],  namely,  the 
angular  frequency,  the  angle  of  incidence  and  the  angle  of  observation. 

This  is  an  interesting  benchmark  as  a  example  of  cavity  scattering,  being 
an  application  area  of  considerable  interest,  and  because  it  is  well  known  that 
even  minor  changes  in  the  scattering  of  a  cavity  can  lead  to  dramatic  changes 
in  the  scattered  fields  due  to  a  variety  of  electromagnetic  phenomena,e.g.,  quasi 
resonant  behavior. 

2.2.  Invisible  Pacman  Benchmark 

The  second  scatterer  we  consider  is  a  perfectly  conducting  2D  cylinder  with 
a  cut-out  wedge.  The  configuration  is  illustrated  in  Figure  1.  Here,  Ow  denotes 
the  angle  of  the  wedge  and  is  a  critical  parameter  for  this  case.  We  will  also  allow 
the  observation  angle  Or  and  angular  frequency  uj  to  vary.  A  curvilinear  PML  [7] 
is  applied  around  P far  to  simulate  the  unbounded  nature  of  the  computational 
domain. 

An  interesting  phenomenon  of  this  problem  is  that  the  scattered  fields  and 
thus  the  RCS  change  dramatically  with  only  a  small  change  in  the  wedge  angle 
Ow  In  Figure  2,  we  show  the  RCS{uj  =  IOtt,  =  0,6^^)  in  dB  when  Ow  = 
0°,18.5°  and  21.5°.  We  note  in  particular  the  22dB  drop  in  the  monostatic 
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Figure  1:  Geometry  of  the  computational  domain  for  the  open  cavity  problem  (left)  and 
invisible  pacman  problem  (right). 

scattering  (when  0^  =  0)  with  just  3  degrees  of  change  in  Ow  rendering  the 
Pacman  substantially  less  visible.  One  of  the  challenges  in  this  work  is  to 
demonstrate  that  reduced  basis  methods  are  capable  of  capturing  such  critical 
behavior  across  variations  in  geometry  and  frequency. 


Figure  2:  Radar  cross  sections  for  the  Pacman,  lOlogio(RCS),  versus  observation  angle  Or- 
Three  cases  with  different  wedge  angles  6w  are  plotted.  In  the  three  cases  6i  =  0  and  cn  =  IOtt. 


2.2.1.  Geometric  Transformation 

To  apply  the  RBM  with  computational  advantage,  the  parametric  problem 
needs  to  be  written  on  a  fixed  domain.  To  address  this,  we  transform  the 
problem  on  a  domain  depending  on  the  angle  Ow  to  a  parametric  problem  on 
a  fixed  domain  through  a  map  from  a  reference  domain  with  a  fixed  0\y  to  a 
physical  domain  with  a  general  Ow 

Let  (x,  y)  and  (p,  0)  be  rectangular  and  polar  coordinates  on  the  reference 
domain,  and  (x,  y)  and  (p,  0)  be  those  on  the  physical  domain.  The  mapping  of 
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the  domain  is  achieved  through  a  sequence  of  mappings  with  the  first  and  last 
being  the  standard  transformation  between  polar  and  rectangular  coordinates: 


X 

y 


P  =  P 

0  =  fm{0) 


(9) 


The  function  fm{')  in  the  second  mapping  is  a  continuous,  piece- wise  smooth 
function  that  is  one-to-one  and  such  that  /m([~7r,  —O^y])  =  [— tt,  —Ow]i  ^w]) 

[-0w,0w]  and 

Let  us  first  discuss  the  systems  of  equations  on  the  physical  domain.  We  omit 
the  T  on  all  the  variables  and  operators  for  simplification  of  notations.  Under 
the  assumption  of  /i  =  1,  e  =  1,  the  TM-polarized  system  of  time  harmonic 
Maxwell’s  equations  (5)  can  be  rewritten  using  polar  coordinate  as 


iujEz 
\  iujHp 
iuoHQ 

\ 


_  1  (d{pHe)  _  dH,\ 
p  y  dp  89  J  ^ 

—  1  dE, 

p  89  ’ 

_  8^ 

8p  • 


(10) 


Next,  we  incorporate  the  perfectly  matched  layer  to  truncate  the  computational 
domain.  To  achieve  this,  we  follow  [7]  for  the  derivation  under  curvilinear 
coordinate  as  a  change-of- variable  technique.  For  this,  we  define 


P  = 


p>a, 
p  <  a. 


Applying  the  PML  layer  is  then  equivalent  to  replacing  all  p  in  (10)  by  p.  Here, 
cr  is  a  piecewise  quadratic  function  of  p  (constant  along  the  6>— direction). 
It  is  identically  zero  in  the  non-PML  region  and  monotonically  increasing  along 
the  p— direction  from  the  PML/non-PML  interface  to  the  exterior  boundary. 
To  simplify  the  system,  we  let 

p>a, 
p  <  a, 

and  obtain,  after  some  manipulation,  the  polar  coordinate  system 

<  {iuj  +  a)Hp  =-l^,  (11) 

{iui  +  a)He  =%-■ 

Converting  to  rectangular  coordinate  we  recover,  after  some  algebraic  manipu¬ 
lations. 


—  {iuj  T  cr )  {iuj  -|-  (7 ^E2;  T  AV  • 


Hy 


=  0, 

=  0. 


(12) 
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Here, 


(a  —  a)  sin6>cos6> 
icj  -j-  O'  cos^  0  a  sin^  0 


=  (CT  -  cr) 


■iu)  —  a  sin^  9  —  a  cos^  9 
{a  —  a)  sin6>cos6> 

sin  0  cos  0  sin^  0  \  f  0 

0  —  sin  cos  0  /  ^  I  iu  ^  a 


—  cos 


-tuj  —  a 
0 


Finally,  we  obtain  the  parametric  system  on  the  reference  domain  to  which 
we  apply  the  reduced  basis  method 


—  (^iuj  “h  (7 ) (^iu)  “h  (J “h  A.  J  ^ V  ' 


ir-T 


VE, 


Here,  JJ  is  the  Jacobian  matrix  of  the  mapping. 


H,. 


=  0, 

=  0. 


(13) 


given  as 


J  = 


d{x,  y)  ^  1 

d{x,  y) 


XX  +  finVy  xy  -  f^xy 

xy  -  f'mxy  yy  +  f'mxx 


3.  The  reduced  basis  method 

In  the  following,  we  outline  the  procedure  applied  to  build  the  reduced  basis. 
3.1.  The  truth  approximation 

We  first  describe  the  finite  element  method  to  solve  the  parametric  systems 
(5)  and  (13)  that  we  want  to  approximate  for  a  particular  parameter  choice. 
We  let 


and,  for  a  given  a  mesh  7^,  we  use  a  discontinuous  Galerkin  method  [11]  to 
solve  both  equations  as  in  [6]  but  without  the  elimination  of  E^  since  this  can 
not  be  done  when  the  first  order  absorbing  boundary  condition  or  the  PML 
are  present.  Here,  we  omit  the  details  (see  e.g.  [6,  11]  for  the  details)  of  the 
formulation  and  simply  denote  the  resulting  bilinear  form  as 

The  truth  approximation  for  the  scattered  field  u-^{u)  is  then  the  solution  of 
a-^(u-^(6>w),  v;^w)  =  v;^w),  V  v  G  X-^ .  (14) 
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Here,  we  define  the  following  finite  element  space 

=  {v  G  {L‘^{Th)f  :  for  all  elements  K  G  Th,v\K  G  {Pk{K)f},  (15) 

and  is  the  incident  field  —d2  . 

Remark  1.  The  volume  terms  on  the  right  hand  side  of  (14)  involves  the  in- 
eident  field  and  ean  he  removed  in  praetiee  sinee  it  automatieally  satisfies  the 
three  first  equations  in  (5)  at  the  eontinuous  level. 

3.2.  Construetion  of  the  redueed  basis.  The  greedy  algorithm 

The  next  step  in  the  approximation  is  to  replace  the  approximation  space 
(15)  for  the  DG  formulation  (14)  by  a  lower  dimensional  space  of  the  form  (see 

(3)) 

=  span  {u-^{iyi),  i  G  {1, . . . ,  A^}}  .  (16) 

As  indicated  previously,  the  selection  of  these  parameters  z/^,  i  G  {1, . . . ,  A^}  is 
crucial  to  achieve  good  approximation  properties  of  the  reduced  space.  This 
selection  is  performed  through  a  greedy  algorithm,  presented  for  completeness 
in  Figure  3  (see  [17,  20,  6]  for  the  details).  The  key  point  of  the  algorithm  is 
line  3  where  the  truth  approximation,  i.e.,  solution  of  the  problem  (14)  or  (2), 
is  compared  with  the  reduced  basis  solution,  recovered  through  (3)).  The  value 
of  the  parameter  for  which  a  bound  for  this  error  is  maximized  is  taken  to  be 
the  next  parameter  values. 

In  practice,  the  set  V  is  replaced  by  a  discrete  training  set  .  Moreover, 

the  a  posteriori  error  estimator  evaluated  at  low  computational  cost  , 

is  assumed  to  be  a  rigorous  and  optimal  upper-bound  of  the  error  [17,  20,  6]. 


1 

2 

3 

4 

5 


=  span  {iz-^(z/i)}  (with  ui  arbitrarily  chosen)  ; 
for  A^  =  1, . .  .,Nmax  do 
Z/7V  =  argmax  A7v(z^); 

uED 

=  span  z  G  {1, . . . ,  A/'}}; 

end 


Figure  3:  Algorithm  to  build  the  reduced  basis  space 
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3.3.  The  empirical  interpolation  method 

The  bilinear  and  linear  forms  in  our  benchmark  problems  both  contain  func¬ 
tions  that  are  non-affine  with  respect  to  the  parameter.  For  an  example,  the 
linear  form  involved  in  the  computation  of  the  output  (7)  contains  terms  of  the 
form 

j,  i>{x)  g{x-,v)dC,  (17) 

where 

g{x^  v)  =  exp(-ico’/c0^  •  x) 

=  co^{uj{xco^0r  ^  y^vnOr))  —  i  sin(co’(x  cos6>^  +  ^sin6>^)) . 

" - V - "  " - V - " 

gi{x,v)  g2{x,v) 

(18) 

The  affine  assumption  enables  the  offline-online  computational  strategy  and  is 
hence  a  key  part  of  computational  efficiency  of  the  scheme.  When  violated,  one 
idea  is  to  approximate  the  function  *)  in  (17)  by  a  linear  combination  of 
functions  with  separate  dependence  of  the  parameter  variable  and  the  spatial 
variable,  i.e., 

M 

g{x;  v)  K,  gM{x]  t/)  =  ^  (09m(®), 

m=l 

using  the  empirical  interpolation  method  [2,  10].  This  leads  to  an  approximate 
linear  form  clearly  satisfying  the  affine  assumption.  The  empirical  interpolation 
method  is  based  on  i)  the  expansion  of  the  function  •)  over  a  space  accounting 
for  the  parameter  dependence  and  ii)  a  robust  interpolation  procedure.  For 
completeness  of  the  paper  we  briefly  present  these  two  steps  [2,  10]. 

Construction  of  the  approximation  space.  To  build  an  approximation  space  well 
adapted  to  the  parameter  dependence  we  consider 

:=  span  {g{--,  1  <  m  <  M) . 

The  set  of  parameters  =  {z/^,  1  <  m  <  M}  is  built  hierarchically  using  the 
greedy  algorithm  in  Figure  4. 

Given  a  set  of  parameters  (and  its  associated  approximation  space  VF^), 
the  next  parameter  is  selected  as  the  one  where  the  difference  between 

and  its  best  approximation,  measured  through  a  suitable  norm  ||  •  ||, 

5m(-,m)  :=  argniin  \\g{-,  fi)  -  z\\  ,  (19) 

is  maximized  in  In  practice,  the  set  V  on  line  4  is  replaced  by  a  training  set 

'^train  ^  addition,  we  use  a  L^-type  norm  in  such  a  way  that  computing 

(19)  amounts  to  the  solution  of  a  linear  system. 
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1  =  0; 

2  =  {0}; 

3  for  M  =  1, . . . ,  Mmax  do 

®  = '^M-l 

6  =  span  ig{-,  v),  v  G 

7  end 

Figure  4:  Algorithm  to  build  the  approximation  space 


1  n  =  0; 

2  for  M  =  1, . . . ,  M^ax  do 

3  Compute  gM-i{-\v\^)  G  W^m-i  interpolating 

5(-;  ^m)  oil  the  points 

4  Compute  rM(-)  =  9{-]Vm)  - 

=  arg  esssup|rM(-)l,  =  i^M-i  FmI; 

^  tcGF 

6  9m(-)  =  rM{-)/rM{xlj); 

7  end 

Figure  5:  Algorithm  to  compute  the  interpolation  points 
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Setting  up  the  interpolation  proeedure.  The  interpolation  points 

l<m<M}, 

that  are  selected  with  the  algorithm  in  Figure  5  lead  to  an  efficient  and  ro¬ 
bust  interpolation.  In  addition,  this  procedure  provides  a  new  basis  ^  ^ 

{!,•••  ,  M}  for  the  interpolation  space  such  that  the  computation  of  the  inter¬ 
polation  coefficients  (u)  amounts  to  solving  a  triangular  linear  system. 

Notice  that  the  family  of  interpolation  points  is  hierarchical.  Moreover,  the 
algorithm  is  well  defined  unless  tmC*)  7^  0.  This  would,  however,  imply  that  the 
interpolation  procedure  would  be  exact  for  all  n  e  V  when  using  M  —  1  terms 
on  the  expansion. 

In  practice,  the  search  of  the  interpolation  point  in  line  5  is  not  performed 
over  r  but  on  a  finite  number  of  points  included  in  T  (usually  the  nodes  of  a 
mesh).  Moreover,  the  two  algorithms  in  Figures  4  and  5  are  interwined,  see 
[2,  10].  That  is,  one  more  basis  function  gM(')  is  added  in  each  step  together 
with  the  next  interpolation  point  Lastly,  the  best  approximation  error  in 
(19)  is  replaced  by  interpolation  error  as  shown  by  line  3  of  Figure  5. 


4.  Numerical  Results 

In  this  section,  we  discuss  the  numerical  results  for  the  two  benchmark  prob¬ 
lems,  seeking  to  confirm  the  high  efficiency  and  good  accuracy  of  the  reduced 
basis  method  applied  to  complex  scattering  problems. 

4.1.  Open  Cavity  Benehmark 

The  governing  equations  are  solved  using  the  discontinuous  Galerkin  method 
with  third  order  polynomials  over  a  nonconforming  mesh  locally  refined  toward 
the  corners  and  tips  of  the  scatterer.  The  total  number  of  unknowns  is  7520.  In 
Figure  6,  we  illustrate  the  fields  when  cj  =  27r  and  Oi  =  7r/4. 

The  linear  functional  associated  with  the  output  of  interest  does  not  satisfy 


Figure  6:  Solution  of  the  problem  on  the  open  cavity  (from  left  to  right:  Ez,  Hx  and  Hy)  for 
uj  =  27t  and  6  =  A. 
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the  affine  assumption  as  discussed  in  Section  3.3).  For  this  reason  we  begin  by 
checking  the  performance  of  the  EIM  used  to  interpolate  the  function  gi{'r)  in 
(18)  (the  results  for  g2{'r)  are  similar)  to  recover  an  affine  linear  form. 

4.1.1.  Empirical  Interpolation  Method 

The  parametric  space  is  in  this  case  two  dimensional  with  v  =  (cu,dr)  since 
di  only  appears  in  the  right-hand  side  of  the  PDF.  Thus  we  need  to  interpolate 
gi(’;  ’)  on  the  boundary  T  where  the  output  is  computed.  The  approximation 
space  for  this  interpolation  is  obtained  using  the  I/^(r)-norm.  Hence,  the  com¬ 
putation  of  the  best  approximation  (see  (19))  amounts  to  solving  a  linear  system 
with  size  equal  to  the  number  of  degrees  of  freedom  on  the  boundary  T.  The 
interpolation  points  are  computed  using  the  algorithm  in  Figure  5,  where  the 
maximum  is  sought  (for  simplicity)  among  the  nodes  used  by  the  approximation 
and  not  in  the  entire  F. 

1-B  parameter  experiment.  To  begin,  we  only  consider  one  parameter,  the  ob¬ 
servation  angle  Or  ^  V  :=  [0,27r],  while  keeping  the  angular  frequency  fixed. 
The  empirical  interpolation  method  amounts  in  this  case  to  the  approximation 
of  a  plane  wave  with  angular  frequency  uj  and  an  arbitrary  measurement  an¬ 
gle  Or  using  a  finite  number  of  plane  waves  with  the  same  angular  frequency. 
The  results  are  displayed  in  Figure  7  for  cj  =  O.Stt  and  uj  =  2.57r.  The  greedy 
search  performed  on  the  line  4  of  the  algorithm  in  Figure  4  is  done  over  the 
training  set  Strain  =  {27r/c/360}|^Q.  The  maximum  dimension  of  the  dis¬ 
crete  space  Mmax  is  the  smallest  integer  such  that  the  best  approximation  error 
^m(^)  •=  is  less  than  e  =  2.5  x  10“^"^  uniformly  on  Strain- 

In  the  first  column  of  Figure  7  we  illustrate  the  parameters  that  are  se¬ 

lected  by  the  algorithm,  hence  defining  the  approximation  space.  In  the  second 
column  we  plot  the  corresponding  interpolation  points  that  are  selected 

by  the  algorithm  in  Figure  5  for  both  configurations. 

To  test  the  quality  of  the  approximation  space  and  the  greedy  interpo¬ 
lation  procedure,  we  consider  a  random  set  of  1000  parameters  Smci  uniformly 
distributed  in  V.  We  compute 

max  and  niax  ||5(-;  0  -  (20) 

ve^rnd 

and  show  these  in  the  third  column  of  Figure  7.  Note  that  the  first  quantity 
is  the  maximum  error  committed  by  the  best  approximation  (projection)  over 
E^rnd  while  the  second  quantity  is  the  maximum  error  in  which  case  the  inter¬ 
polation  is  used  over  'Emd-  Both  quantities  show  exponential  decay  when  the 
number  of  elements  is  increased.  Moreover,  the  quality  of  the  projection  and 
the  interpolation  is  comparable,  implying  that  the  interpolation  points  are  well 
selected  and  the  interpolation  procedure  is  robust.  As  expected,  for  a  higher 
angular  frequency  we  need  more  terms  in  the  approximation  to  achieve  a  similar 
accuracy. 
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Int.  points  for  cos(a)(x  cos  0  +  y  sin  0)).  co  =  0.5  ti 


0)  =  0.5  7t 


Figure  7:  Empirical  interpolation  of  the  function  gi{x]  =  cos(u;(x cos(^r)  +  ysinOr)) 

with  6r  G  [0,  27t]  when  uu  =  O.Stt  (top)  and  u  =  2.507r  (bottom)  on  F.  In  the  first  (resp.  second) 
column,  we  represent  the  set  (resp.  the  interpolation  points  ).  The  larger  the 

point  is  the  earlier  it  has  been  chosen.  Plotted  in  the  third  column  are  max^y^s^^^  Il5'(s  “ 
5'm(-;^)||  and  max^es^^d  Il5'(-;^)  -5'm(-;^)||- 


2-D  parameter  experiment.  We  interpolate  the  function  for  the  case 

where  both  {uj,0r)  varies  on  the  set  V  =  [0,  Stt]  x  [0,27r].  In  this  case  the 
empirical  interpolation  method  approaches  a  plane  wave  with  arbitrary  {uj^Or) 
in  the  given  range  by  a  linear  combination  of  a  finite  number  of  plane  waves. 
The  results  are  displayed  in  Figure  8.  On  the  left  we  represent  the  parameters 
selected  by  the  greedy  process  on  the  line  4  of  the  algorithm  in  Figure  4,  run  over 
the  training  set  Strain  =  {27r/c/360}|^Q  x  {37r/c/100}^^Q.  Note  that  the  high 
frequencies  (where  the  function  oscillates  more)  are  preferred.  In  the  middle  we 
represent  the  interpolation  points  selected  by  the  algorithm  in  Figure  5. 

To  test  the  quality  of  the  approximation  space  and  the  interpolation  points 
we  evaluate  the  algorithm  over  a  set  of  40000  parameters  randomly  chosen  in  the 
set  [0,27r]  x  [0,37r],  and  compute  the  quantities  in  (20).  The  results  are  plotted 
in  Figure  8  on  the  right.  We  observe  that  with  15  terms  the  error  begins  to 
decay  exponentially.  With  about  25  terms  a  10“^  accuracy  is  achieved.  We  also 
observe  that  the  larger  the  range  of  frequencies  is,  the  higher  the  dimension  of 
the  approximation  space  is  needed.  Hence,  in  applications  where  large  ranges 
are  required,  a  computational  strategy  for  splitting  this  range  into  smaller  ones 
to  control  overall  online  cost  is  interesting  and  perhaps  even  essential. 

In  the  next  section,  when  coupling  the  EIM  with  the  RBM,  we  use  a  large 
number  of  terms  on  the  expansion  (M  =  35)  in  order  to  ensure  a  good  approx¬ 
imation  uniformly  across  the  parameter  space. 

4.1.2.  RBM  with  EIM 

Since  the  non-affinity  only  resides  in  the  linear  forms  for  this  problem,  we 
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Selected  angles  for  cos(a)(x  cos  0  +  y  sin  0)) 
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Interp.  points  for  cos(a)(x  cos  0  +  y  sin  0)) 


Figure  8:  Empirical  interpolation  of  the  function  (uj,6r))  =  cos  {u(x  cos  (Or)  +  ysin^^’)) 

with  (cn,  Or)  G  [0,  Stt]  x  [0,  27r]  on  F.  In  the  first  (resp.  second)  column  we  represent  the 
set  (resp.  the  interpolation  points  ).  The  larger  the  point  is  the  earlier 

it  has  been  chosen.  Plotted  in  the  third  column  are  \\g(-]u)  —  cind 

max^es^^d  Il5'(-;^)  -  9m{-;i^)\\- 

only  need  to  apply  the  empirical  interpolation  method  on  the  linear  forms. 
Once  we  have  obtained  an  affine  expansion  of  the  linear  forms,  the  next  step 
of  the  offline  computations  consists  of  computing  the  reduced  basis  that  will 
be  used  in  the  online  part.  We  recall  that  the  output  of  interest  depends  on 

three  parameters  (uu^Oi^Or)  G  P  =  [7r/2,57r/2]  x  [0,  27r]  x  [0,27r].  We  start  by 

splitting  the  parameter  domain  [8]  into  three  sub-domains  (corresponding  to 
low  frequencies,  medium  frequencies  and  high  frequencies)  V  =  where 


TT  Sir 

Sit 

57r 

■Di  = 

2’  Y 

x[0,27r]^,  I>2  = 

x[0,27r]^  and  P3  = 

27r,-^ 

This  allows  us  to  use  different  discretization  spaces  for  the  truth  approximation 
(for  small  values  of  cj  a  coarser  mesh  could  be  used)  depending  on  the  parameter 
sub-domain.  We  apply  the  greedy  algorithm  3  over  the  train  sets 


'Strain 

r  2l7r 

1  25  fl 

'Strain 

^2 

95 

r  Stt  riTT  1  ^ 

~  12“^  2-25  /n=0 

'Strain 

^3 

=  {2»+»)r.o 

x{|f 

}?Io  X  {SP}L.  C 

Note  that  the  problem  is  more  sensitive  to  parameter  changes  for  large  values  uj. 
The  density  of  the  train  set  is  selected  accordingly.  In  addition,  the  construction 
of  the  reduced  basis  is  independent  on  which  appears  only  in  the  linear  form 
associated  to  the  output.  This  allows  us  to  run  the  greedy  algorithm  over  two- 
dimensional  sets. 

In  Figure  9,  we  show  the  parameters  selected  by  the  greedy  algorithm  in 
the  three  cases  (left),  and  the  evolution  of  the  maximum  over  the  train  sets 
^trazn,  ^  ^  {1,2,3}  of  the  a  posteriori  error  estimator  (right).  We  observe 
exponential  convergence  of  the  a  posteriori  error  estimator  and  see,  as  expected, 
that  the  rate  of  convergence  is  lower  when  we  consider  larger  parameter  sets  or 
high  values  for  uj. 

Let  us  now  test  the  quality  of  the  reduced  basis  space  in  the  online  stage. 
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Figure  9:  Construction  of  the  reduced  basis  when  (u;,0i)  is  in  [7r/2,37r/2]  x  [0,  27r]  (in  blue), 
[37r/2,27r]  x  [0,  27r]  (in  red)  and  [27r,57r/2]  x  [0,  27r]  (in  green).  On  the  left  we  represent  the 
parameters  selected  by  RBM  (the  larger  the  markers  the  earlier  they  have  been  chosen).  On 
the  right  is  shown  the  evolution  of  the  maximum  over  the  train  sets  i  G  {1,2,3}  of 

the  a  posteriori  error  estimator  when  the  dimension  of  the  reduced  basis  is  increased. 


MonoStatic  RCS,  Figures  10  and  11.  We  start  by  computing  the  monostatic 
RCS  in  dB  (i.e.,  we  consider  Or  =  Oi)  with  {uj^Oi)  G  [7r/2,57r/2]  x  [0,27r].  The 
results  are  shown  in  Figure  10  (left).  The  logarithm  in  base  10  of  the  a  posteriori 
error  estimator  for  the  linear  output  F(',  •,  •)  is  also  plotted  (right).  Note  that 
depending  on  the  sub-domain  in  which  the  parameter  falls  (Pi,  V2  or  P3),  a  dif¬ 
ferent  reduced  basis  space  is  used.  For  both  plots  we  perform  the  computations 
over  the  following  grid 


'^test 


'^test 

^2 


^test 

^3 


fTT  I  ^  r/2Z7r  2Z7r\'|360 

I  2  ^  100  /  n=0  I  V  360  ’  360  /  /  /=0  ^  ^  ’ 

/  UTT  1  V  /  f  2Z7r  21tt  \  1  ^  'j^ 

I  2  ^  2-100 /n=0  IV360’  360>'/z=0  ^ 


which  is  not  a  subset  of  the  training  set.  We  have  used  190,  220  and  250  elements 
for  the  computations  over  and  respectively.  We  present,  in 

Figure  11,  the  monostatic  RCS  for  uj  G  {1.5708,6.2832,7.2257}  computed  with 
two  different  number  of  elements  on  the  reduced  basis.  We  can  observe  that  in 
both  cases  the  approximation  is  very  good  even  if  the  output  is  less  reliable  with 
a  low  number  of  elements  on  the  basis.  Adding  a  few  elements  to  the  reduced 
basis  space  dramatically  reduces  the  width  of  the  confidence  interval.  Note  that 
these  plots  correspond  to  three  different  cuts  of  the  monostatic  RCS  shown  in 
Figure  10  for  the  values  of  uj  selected  above. 

We  plot,  in  Figure  12,  the  bistatic  RCS  in  dB  for  uj  =  7.4693  and  all  possible 
angles  of  incidence  and  measurement  and  the  logarithm  in  base  10  of  the  a 
posteriori  error  estimator  for  the  output.  For  both  plots  we  have  performed  the 
computations  over  the  following  grid 

s-  =  {7.4693} 
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Figure  10:  Monostatic  RCS  (i.e.  Or  =  Oi)  in  dB  (left)  for  {uo.Oi)  G  Vi  U  V2  U  with 
Vi  =  [7r/2,37r/2]  x  [0,  27r],  V2  =  [37r/2,27r]  x  [0,  27r]  and  =  [27r,57r/2]  x  [0,27r].  For  each 
of  these  three  regions  we  compute  a  different  reduced  basis.  Logarithm  in  base  10  of  the 
estimator  for  the  output  of  interest  is  on  the  right.  The  computations  have  been  done  with 
a  reduced  basis  of  dimension  equal  to  190  for  (u;,0i)  G  Vi,  220  for  (cj,0i)  G  T>2  and  250  for 
(cu,  Oi)  G  'Ds- 


that  is  not  a  subset  of  the  train  set.  In  Figure  13  we  present  the  monostatic 
RCS  and  the  bistatic  RCS  in  dB  for  Oi  G  {0.0000, 1.2566}.  We  observe  that  the 
convergence  of  the  confidence  region  provided  by  the  estimator  rapidly  shrinks 
when  the  dimension  is  increased.  Note  that  these  plots  correspond  to  the  re¬ 
striction  of  the  plot  12  to  the  solid,  dashed  and  dotted  lines. 
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=  1.5708.  Dim  =  150 


0. 


CO  =  1.5708.  Dim  =  180 


0. 


Figure  11:  From  left  to  right,  monostatic  RCS  in  dB  for  uj  G  {1.5708,6.2832,7.2257}.  In 
black,  the  output  computed  with  the  reduced  basis  method.  The  blue  region  is  the  confidence 
region  provided  by  the  a  posteriori  error  estimator.  In  the  first  row,  we  use  150,  150,  and  180 
bases  respectively  and  the  reduced  basis  output  already  provides  a  good  approximation  of  the 
actual  output.  However,  the  confidence  region  for  some  angles  is  rather  wide.  In  the  second 
row,  we  enrich  the  reduced  basis  (from  150  to  180  and  from  180  to  220)  to  certify  the  output. 


M  RCShnJia  Dim  >330 


Figure  12:  Bistatic  RCS  in  dB  (left)  for  {6i,6r)  G  [0,27r]^  and  uo  =  7.4693,  and  logarithm  in 
base  10  of  the  estimator  for  the  output  of  interest  (right).  The  computations  are  done  with  a 
RB  dimension  equal  to  220. 
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0J=  7.4693  Dim  =  180  co  =  7.4693.  9.  =  0.  Dim  =  1 80  co  =  7.4693.  9.  =  1 .2566.  Dim  =  180 


Figure  13:  From  left  to  right,  monostatic  RCS,  bistatic  RCS  in  dB  with  6i  =  0.0000,  and 
bistatic  RCS  with  6i  =  1.2566  in  dB.  u  =  7.4693  in  all  cases.  Plotted  in  black  is  the  output 
computed  with  the  reduced  basis  method  using  180  (top)  and  220  (bottom)  bases.  The  blue 
region  is  the  confidence  region  provided  by  the  a  posteriori  error  estimator.  The  RB  dimension 
is  180  (top)  and  220  (bottom). 
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^.2.  Invisible  Pacman  Benchmark 

For  the  Pacman  scattering  problem,  we  consider  two  parameters:  angular 
frequency  uj  and  the  wedge  angle  Ow^  and  apply  a  reduced  basis  method  to  the 
parametric  equations  (13)  derived  in  Section  2.2.  To  emphasize  the  dependence 
on  the  parameters,  the  radar  cross  section  (RCS),  equations  (7)  and  (8),  is 
denoted  by  RCS{uj^Ow,^i^^r)-  For  an  example,  RC5'(37r,  18.5°,  0,  in  the 
caption  of  a  figure  means  that  the  plot  is  the  RCS  versus  observation  angle  Or 
with  frequency  37r,  wedge  angle  18.5°,  and  incident  angle  0. 

The  scatterer  has  a  fixed  diameter  (5  wavelengths  when  uj  =  IOtt)  and  we 
have  used  a  curvilinear  PML  [7]  of  fixed  thickness,  corresponding  to  2.5  wave¬ 
lengths  away  from  the  scatterer,  in  an  annulus  of  width  1.5  wavelengths  when 
UJ  =  IOtt.  The  PML  has  a  maximum  damping  coefficient  60. 

4-2.1.  Truth  approximation 

We  begin  by  studying  some  aspects  of  the  truth  solver  for  the  problem  (13). 
We  fix  0\y  =  18.5°,  the  reference  angle  for  the  wedge,  throughout  the  paper. 

When  Ow  =  0,  we  have  the  exact  solution  to  compare  our  result  with.  Here, 
we  solve  the  scattering  problem  for  a  cylinder  with  two  frequencies:  Sir,  and  IOtt 
and  calculate  the  RCS.  We  see,  from  Figure  14,  that  the  finite  element  solver  is 
highly  reliable  with  relative  error  on  the  level  of  10“^  for  higher  frequency  and 
10“^  for  lower  frequency. 


Figure  14:  RCS  for  the  cylinder  10logio(RCS((jj,0,0,Or))  with  wave  number  a;  =  Stt  for  the 
left  and  IOtt  for  the  right.  From  top  to  bottom:  exact  solution,  approximation,  relative  error. 


To  test  the  efficiency  of  this  PML  applied  to  this  problem,  we  compare 
the  RCS  with  that  obtained  by  using  only  a  first  order  absorbing  boundary 
condition,  Silver-Miiller  boundary  condition.  It  is  evident  from  Figure  15  that 
the  application  of  the  curvilinear  PML  improves  the  results.  Here,  the  reference 
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30 


Figure  15:  Left:  comparison  of  the  bistatic  RCS  {6w  =  18.5°),  10/o5fio(RCS'(107r,  18.5,  0,  ^r)), 
with  application  of  curvilinear  PML  and  with  a  first  order  absorbing  boundary  condition. 
Right:  relative  error  of  the  computation  with  curvilinear  PML. 


is  a  proven  solver  for  the  Pacman  scatterer  on  a  square  domain  with  rectilinear 
PML  and  no  geometric  mapping. 

Next,  we  study  the  effect  of  three  possible  forms  of  mapping  /m(')  in  (9). 
We  show  here  the  resulting  monostatic  RCS  in  dB  from  the  parametric  solver 
with  Ow  =  21.5.  In  Figure  16,  “P2P1P2”  indicates  that  the  polynomial  de¬ 
grees  of  the  mapping  on  [— tt,  — ^^],  [—0\y^^w]  ^ 

respectively;  similarly  “P2P2”  degrees  on  [— 7r,0]  and  [0,7r];  “Pi  Pi  Pi”  degrees 
on  [-TT, -6*^]  U[6'Jy,7r]  and  [-01^,91^].  Naturally,  we  require  /m(-7r)  =  — tt, 
fm{-9lv)  =  -9w,  fm{9w)  =  9w,  fm{T^)  =  for  Continuity  and  also  demand 
—regularity  if  there  are  additional  degrees  of  freedom  allowing  us  to  do  so 
(first  and  second  cases). 


Figure  16:  Comparison  of  three  different  mappings 


We  observe  that  the  “Pi  Pi  Pi”  mapping  yields  the  best  overall  result.  The 
reason  for  this  is  found  in  the  observation  that  this  map  ensures  the  most  reg¬ 
ularity  in  the  physical  domain.  Hence,  we  shall  use  this  mapping  in  the  truth 
solver. 

4.2.2.  RBM  for  one  parameter  without  EIM 

The  interesting  phenomenon  described  in  Section  2.2  and  exemplified  by 
Figure  2  suggest  the  study  of  the  dependence  of  scattering  on  varying  Ow 
Therefore,  we  fix  cj  =  IOtt,  treat  Ow  as  our  parameter  in  this  section  and 
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develop  a  reduced  basis  method  to  the  Pacman  scattering  problem.  Here,  we 
set  the  parameter  domain  for  Ow  to  be  [8.5°,  28.5°].  Contrary  to  the  problem 
treated  in  Section  4.1  where  only  the  linear  forms  were  non-afhne,  the  bilinear 
form  for  this  problem  is  also  non-affine  thanks  to  its  dependency  on  Ow  through 
the  geometric  transformation.  Note  that  only  the  parameter  Ow  is  responsible 
for  the  non-affinity  since  the  bilinear  form  is  affine  in  uj  (see  equation  (13)). 
EIM  will  be  applied  in  the  next  subsection  where  it  is  shown  that  it  works  well 
for  this  problem  and,  as  expected,  does  not  degrade  the  quality  of  the  RBM 
solution.  However,  the  results  in  this  subsection  are  obtained  without  EIM  for 
simplicity. 

The  procedure  of  building  the  reduced  basis  space  is  standard  and  we  omit 
the  details.  The  accuracy  of  the  reduced  basis  solution  is  certified  by  the 
residual-based  a  posteriori  error  estimate.  The  inf-sup  number  will  not  be¬ 
come  zero  since  this  is  an  open  problem.  The  resulting  error  estimate  is  cheap 
to  obtain  online.  It  also  guides  the  selection  of  the  parameters  and  thus  the 
building  of  the  reduced  basis  space  in  the  greedy  algorithm  in  3.  See  e.g.  [6,  19] 
for  details.  Eor  this  example,  this  algorithm  has  been  run  over  the  train  set 

strain  257. 

In  Eigure  17  we  illustrate  the  29  parameter  instances  that  the  reduced  basis 
method  picks  together  with  their  orders.  The  first  point  is  picked  randomly 
to  start  the  greedy  algorithm.  It  is  also  interesting  to  note  that,  in  order  to 
capture  the  critical  angles,  the  greedy  algorithm  does  not  necessarily  select 
points  clustering  around  those  angles.  This  underscores  the  strength  of  the 
method  and  that  its  ability  to  capture  the  critical  phenomenon  does  not  rely  on 
sampling  around  the  critical  points  in  the  parameter  domain. 

We  plot,  in  Eigure  18,  the  history  of  convergence  of  the  RB  solutions  for  the 
worst  of  the  120  randomly  selected  parameter  values.  Exponential  convergence 
is  clearly  observed.  Moreover,  the  error  estimate  decreases  exponentially  with 
about  the  same  rate.  Thus,  the  effectivity  index  is  roughly  constant  across  the 
many  magnitudes  of  decrease  in  the  RBM  error.  With  only  20  bases,  we  obtain 
an  accuracy  at  the  level  of  10“^.  Using  the  underlying  mesh  of  our  numerical 
experiment,  instead  of  solving  systems  of  dimension  above  200,  000  x  200,  000, 
we  only  need  to  solve  problems  of  dimension  60  x  60. 

Next,  we  study  the  monostatic  scattering  as  a  function  of  the  wedge  angle, 
that  is  10logio{RC S {IOtt ,  Ow^  0,  0))-  In  Eigure  19,  we  plot  the  truth  approxima¬ 
tion  (in  black)  and  the  reduced  basis  approximations  with  up  to  only  11  basis 
elements.  We  see  that  the  curve  with  11  bases  and  the  truth  cuve  are  “roughly 
identical”  to  naked  eyes.  We  also  observe  the  dramatic  effect  of  adding  two 
correct  samples  on  the  output. 

As  described  in  section  2.2  and  shown  here  by  the  monostatic  scattering 
curve  in  Eigure  19,  there  are  several  critical  angles  producing  minimal  mono¬ 
static  scattering.  To  visualize  this  phenomenon,  we  show,  in  Eigure  20,  the 
electric  fields  corresponding  to  Ow  =  14.3°,  18.5°  and  21.5°  degrees.  The  dif¬ 
ference  is  easily  noticeable  around  the  wedge. 

Einally,  the  error  estimate  provides  bounds  for  the  output  computed  by  the 
RBM.  In  Eigure  21,  we  plot  the  RBM  output  together  with  the  error  bounds  for 
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Figure  17:  The  29  the  greedy  algorithm  of  the  RBM  picks  to  build  the  RB  space.  Top: 
the  higher  the  vertical  line,  the  earlier  that  point  was  picked.  Bottom:  the  points  scattered 
on  the  monostatic  scattering  curve  10/o5fio(RCS'(107r,  0,  0)),  the  larger  the  marker,  the 

earlier  it  is  selected. 


Figure  18:  The  worst  case  convergence  history  and  the  corresponding  error  estimate  of  the 
RBM  for  120  randomly  selected  parameter  values. 
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Figure  19:  Using  only  11  bases,  we  can  obtain  a  very  accurate  plot  of  the  monostatic  scattering 
{Oi  =  Or  =  ^)  respect  to  the  wedge  angle  Ow 


Figure  20:  Real  part  (top),  imaginary  part  (middle),  and  module  (bottom)  of  the  electric  field 
with  Ow  =  14.3°  (left),  18.5°  (middle)  and  21.5°  (right). 
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three  different  numbers  of  bases.  We  see  that  RBM  output  with  just  17  bases 
captures  all  essential  behavior  with  a  high  level  of  confidence. 


Figure  21:  Error  bounds  of  the  RB  output,  RCS'(107r,  0,  0),  for  three  different  number  of 

bases. 


4.2.3.  RBM  for  two  parameters  with  EIM 

In  this  section  we  consider  two  parameters  in  the  invisible  Pacman  problem: 
the  angular  frequency  uj  and  the  wedge  angle  Ow  As  pointed  out  in  Section 
4.2.2,  both  the  bilinear  form  and  the  linear  forms  for  this  problem  are  non- affine 
due  to  their  dependency  on  Ow  through  the  geometric  transformation.  The 
non- affinity  attributed  to  Ow  is  handled  by  EIM.  To  do  that,  every  function 
depending  on  Ow  in  a  non-afhne  way  (for  example  all  the  components  of  the 
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Jacobian  matrix  associated  to  the  transformation)  are  expanded  using  a  one¬ 
dimensional  empirical  interpolation  method  [10].  EIM  is  applied  in  a  similar 
situation  in  [18,  13].  The  number  of  expansion  terms  M  is  taken  in  {1,2,3}. 
These  numbers  are  quite  small  but  note  that  it  is  applied  to  one-dimensional 
functions  and  there  are  more  than  60  such  functions  in  the  linear  and  bilinear 
forms.  This  results  in  Qa  and  Qm  around  200.  In  this  case,  a  I/^(0)-norm  best 
approximation  is  considered.  It  is  important  to  note  that  the  RBM  procedure 
still  needs  to  be  applied  to  a  two  dimensional  parameter  domain. 

(a) 


Figure  22:  RBM  with  EIM  for  uu  G  [Stt,  Stt]  and  6w  ^  [8.5°,  28.5°]:  (a).  The  parameter 
values  selected  by  RBM  to  build  the  RB  basis.  The  larger  the  marker,  the  earlier  it  has  been 
selected,  (b).  Convergence  of  the  error  with  respect  to  the  dimension  of  the  RB  space  for 
different  number  of  magic  points,  (c).  Convergence  of  the  error  estimate  with  respect  to  the 
dimension  of  the  RB  space  for  different  number  of  magic  points. 


We  first  test  the  method  for  (uj,0w)  ^  [37r,  Stt]  x  [8.5°,  28.5°].  In  Figure  22 
we  show  the  parameter  points  selected  by  the  RBM,  the  worst-case  convergence 
of  the  approximate  solution  (over  the  train  set)  of  the  RBM  coupled  with  EIM 
and  its  estimate.  Exponential  convergence  is  observed.  Note  that  the  green 
curve  corresponding  to  the  case  M  =  1  is  cut  at  the  starting  point  of  the 
plateau  where  the  interpolation  error  begins  to  dominate.  Next,  we  compute 
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Hmdgm  mnglm 


Hedge  eagle 


Figure  23:  RCS{uj^9w^^:^)  in  dB  computed  by  the  RBM  with  EIM  for  oo  G  [tt,  Stt]  and 
Ow  £  [8.5°,  28.5°].  Top:  3D-plot;  middle:  top  view;  bottom:  point-wise  relative  error  of 
RCS{(jj,  0\Y ,  0,  0). 
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the  monostatic  RCS,  that  is  RCS{uj,  Ow^  O5  0)-  Plotted  in  Figure  23  (top)  is  the 
approximation  given  by  RBM  for  {uj^Ow)  ^  [tt,  Stt]  x  [8.5°,  28.5°].  The  RBM 
approximation  is  obtained  by  using  30  bases  for  [tt,  37r]  x  [8.5°,  28.5°]  and  30  bases 
for  [37r,  57r]  x  [8.5°,  28.5°].  We  plot,  in  the  middle,  the  top  view  of  RCS.  The  two 
regions  in  this  parameter  domain  for  the  critical  configurations  that  produces 
very  small  RCS  is  clearly  identified.  At  the  bottom  is  the  point-wise  relative 
error  between  this  RB  approximation  and  the  truth  approximation.  Notice  that 
the  minimum  contour  level  is  set  to  be  0.01  to  show  that  the  approximation 
is  very  accurate  except  for  small  neighborhoods  of  parametric  configurations 
rendering  the  Pacman  invisible.  In  these  small  neighborhoods,  the  errors  are 
substantial.  We  thus  perform  a  new  set  of  computations  using  40  RB  instead  of 
30.  The  result  is  in  Figure  24  showing  that  the  error  is  indeed  getting  smaller 
and  closely  clustering  around  a  few  points. 

The  results  in  Figures  22  and  23  are  obtained  by  using  second  order  accurate 
discontinuous  Galerkin  method  for  the  efficiency  of  our  calculation.  This  suffices 
for  frequency  in  [tt,  57r]  with  the  resulting  truth  approximations  being  accurate 
enough.  However,  when  the  frequency  is  increased,  we  need  high-order  DC 
scheme  to  produce  more  accurate  truth  approximation.  This  is  clearly  shown 
by  Figure  25  where  the  RB  result  is  converging  perfectly  to  (an  inaccurate)  DG 
result  if  second  order  DG  is  used  for  uj  =  IOtt. 


Figure  24:  Point-wise  relative  error  RCSiuo^Ow in  dB  computed  by  the  RBM  using  40 
bases  for  uj  G  [tt,  Stt]  and  6w  ^  [8.5°,  28.5°]. 


Finally,  we  use  a  fourth  order  DG  and  the  resulting  RBM  (coupled  with 
EIM)  to  consider  the  problem  with  {uj^Ow)  ^  [9.87r,  IOtt]  x  [8.5°,  28.5°].  The 
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RCS  in  dB  RCS  in  dB 


2nd  order  DG 


Figure  25:  Comparison  of  the  accuracy  of  RCS(107t,  6w:  0,  0)  in  dB  with  truth  approximations 
computed  by  a  second  and  a  fourth  order  DG  method  for  6\y  G  [8.5°,  28.5°]. 


Figure  26:  RCS{uj,9w:^:^)  in  dB  computed  by  the  RBM  (with  22  bases)  with  EIM  for 
ca  G  [9.87r,  IOtt]  and  Ow  e  [8.5°,  28.5°]. 
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approximate  RCS{uj,0w,^^^)  is  plotted  in  Figure  26  with  22  bases  used.  A 
slice  of  this  picture  with  uj  =  IOtt  shows  comparable  accuracy  with  the  “RBM 
with  11  Bases”  curve  in  Figure  19.  Although  the  complexity  of  the  geometric 
mapping  and  the  system  being  in  a  vector  form  do  lead  to  a  large  number  of 
terms  for  the  affine  decomposition,  the  resulting  RBM  coupled  with  EIM  does 
provide  highly  accurate  approximations  that  converge  exponentially  to  the  truth 
approximations. 


5.  Conclusion 

In  this  paper,  we  have  developed  and  applied  reduced  basis  method  to  two 
non-trivial  electromagnetic  scattering  problems  for  which  the  radar  signature  is 
an  important  output  of  interest.  As  we  have  discussed  in  detail,  this  develop¬ 
ment  has  required  the  development  and  application  of  a  variety  of  techniques 
such  as  empirical  interpolation,  non-trivial  geometric  mapping,  the  inclusion  of 
perfectly  matched  layers  into  the  formulation. 

We  have  demonstrated  that  in  spite  of  the  need  to  combine  these  different 
elements  to  ensure  computational  efficiency,  exponential  convergence  of  the  re¬ 
duced  basis  solution  toward  the  truth  finite  element  approximation  is  obtained 
across  the  entire  parameter  domains.  We  also  demonstrate  that  even  with  very 
narrow  critical  behavior,  the  greedy  approach  captures  this  behavior  very  well, 
even  without  requiring  a  dense  sampling  close  to  the  critical  parameter  values. 

Even  with  these  advances,  several  challenges  still  he  ahead  to  further  mature 
the  reduced  basis  methods  for  complex  scattering  applications.  These  are  mainly 
related  to  increasing  the  online  efficiency,  seeking  to  demonstrate  substantial 
speedup  with  as  much  as  two  to  three  orders  of  magnitude  as  seen  in  other 
applications.  The  essential  components  to  achieve  this  found  in  reducing  the 
complexity  of  the  affine  expansion  recovered  by  the  empirical  interpolations 
for  geometrical  parametric  variations  and  the  general  computational  challenges 
associated  with  large  parameter  ranges  or  high-dimensional  parameter  domains. 
We  hope  to  report  on  such  developments  in  future  work. 
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